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ABSTRACT 

We have performed a series of simulations of clusters of galaxies on the basis 
of the smoothed particle hydrodynamics technique in a spatially-flat cold dark 
matter universe with f2 = 0.3, A = 0.7, and Hq = 70km/s/Mpc as one of the most 
successful representative cosmological scenarios. In particular, we focus on the 
Sunyaev-Zel'dovich effect in submm and mm bands, and estimate the reliability of the 
estimates of the global Hubble constant Hq and the peculiar velocity of clusters v r . 
Our simulations indicate that fractional uncertainties of the estimates of Hq amount 
to ~ 20% mainly due to the departure from the isothermal and spherical gas density 
distribution. We find a systematic underestimate bias of Hq by ~ 20% for clusters 
z ~ 1, but not at z ~ 0. The gas temperature drop in the central regions of our 
simulated clusters leads to the underestimate bias of v r by ~ 5% at z « and by 
~ 15% at z ~ 1 in addition to the statistical errors of the comparable amount due to 
the non-spherical gas profile. 

Subject headings: Cosmology - Dark matter - Hydrodynamics - Galaxies: clusters of 
- Galaxies: X-rays - Numerical methods 
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1. Introduction 

Clusters of galaxies have been extensively observed in radio, optical and X-ray bands. 
Furthermore recent and future observational facilities in mm and submm bands, such as the 
SCUBA (Submillimeter Common-User Bolometer Array), the Japanese LMSA (Large Millimeter 
and Submillimeter Array) project and the European PLANCK mission are expected to open the 
submm window to observe clusters of galaxies via the Sunyaev-Zel'dovich (SZ) effect (Sunyaev, 
Zel'dovich 1972) in addition to the Rayleigh-Jeans region of the spectrum of the cosmic microwave 
background (CMB) where the SZ temperature decrement is reported for about a dozen of clusters 
(e.g., Rephaeli 1995; Kobayashi, Sasaki, Suto 1996). Since the intensity of the SZ effect does not 
suffer from the (1 + z)~ 4 diminishing factor unlike X-ray surface brightness, observations in mm 
and submm bands are much more advantageous for clusters, especially at high z, than those in 
optical and X-ray bands (Barbosa et al. 1996; Silverberg et al. 1997; Kitayama, Sasaki, Suto 
1998). 

As extensively discussed in previous literatures (Silk, White 1978; Sunyaev, Zel'dovich 
1980; Birkinshaw, Hughes, Arnaud 1991; Rephaeli, Lahav 1991; Inagaki, Suginohara, Suto 1995; 
Kobayashi, Sasaki, Suto 1996; Holzapfel et al. 1997; Kitayama, Sasaki, Suto 1998), one can then 
combine these multi-band observations of high-z clusters to determine the cosmological parameters 
and the peculiar velocity of clusters. These procedures, however, usually assume that the gas of 
clusters is isothermal and spherical, while all the observed clusters do exhibit a departure from the 
assumption to some extent. The departure would hamper the reliable estimates of, for instance, 
the Hubble constant H$ and the peculiar velocity v r . To address the question quantitatively, 
we carry out a series of numerical simulations of clusters. We extract simulated clusters both 
at z ps 0.0 and z ps 1.0, and perform the "simulated" observations in X-ray, mm and submm 
bands. Finally we combine the multi-band information to evaluate the statistical and possible 
systematic errors of the estimates of Hq and v r . Our present work extends the previous studies 
of this methodology (Inagaki, Suginohara, Suto 1995; Roettiger, Stone, Mushotzky 1997) for Hq 
and examines uncertainties of the peculiar velocity field as well, paying attention to the projection 
effect and the evolution of clusters. 



2. Numerical Simulation 

As a representative cosmological model, we consider a cold dark matter cosmogony with 
Uq = 0.3, A = 0.7, h = 0.7, cr 8 = 1-0, n = 1 and Q b = 0.015/i~ 2 , where Qq is the density 
parameter, Ao is the cosmological constant, h is Hq i n units of lOOkm/s/Mpc, as is the amplitude 
of the density fluctuation, n is the spectral index of the primordial fluctuations, and 0,^ is the 
baryon density parameter. This model satisfies both the COBE normalization and the cluster 
abundance (White, Efstathiou, Frenk 1993; Kitayama, Suto 1997), and also is consistent with the 
recent indications of no significant evolution of galaxy clusters (see Kitayama, Sasaki, Suto 1998 
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for details). 



The initial conditions of our simulations are generated using the COSMICS package by E. 
Bertschinger and the resulting gas and dark matter particle distributions are evolved with the 
publicly available AP 3 M-SPH code (Hydra) by Couchman et al. (1995). The effects of radiative 
cooling and heating are neglected. We use the ideal gas equation of state with an adiabatic index 
7 = 5/3. Gravitational forces are softened with a physical softening length of 6q = 39/i _1 kpc. 

Our simulations proceed in two steps; first we carry out three low-resolution simulations 
with iV = 64 3 particles each for gas and dark matter in a comoving periodic cubic of L^ ox = 
100/i _1 Mpc (one realization) and 200/i~ 1 Mpc (two realizations). Then we identify clusters at 
z = using the friend-of-friend algorithm with a bonding length of 0.2 times the mean particle 
separation. The virial mass of clusters is computed from all the particles within the virial radius, 
r v i r , from the center of each cluster. The latter is defined so that the mean density inside becomes 
pvir — 187r 2 r2Q 4 /5 c (z), the virialized density at z predicted in the spherical non-linear model 
(e.g., White, Efstathiou, Frenk 1993; Kitayama, Suto 1996), where p c (z) is the critical density 
of the universe at redshift of z. The resulting mass function for the identified clusters from the 
low-resolution simulations is consistent, within a factor of two, with the theoretical Press-Schechter 
mass function (Press, Schechter 1974). Considering a somewhat simplified identification scheme 
adopted here, the agreement is satisfactory and implies that our low-resolution simulations provide 
fairly homogeneous and unbiased catalogues of clusters. 

From the cluster catalogues constructed by the three low-resolution simulations, we select 
nine clusters with mass greater than 1O 14 M0. Then we set a box of side 25/i~ 1 Mpc (5 clusters), 
50/i -1 Mpc (3 clusters) or lOOft^Mpc (1 cluster) at the center of each cluster and fill the box 
with N = 64 3 particles each for gas and dark matter which are assigned the same amplitude and 
phases of perturbation waves as in the initial condition of the low-resolution simulation at z = 32. 
The nine initial conditions for high-resolution simulations are evolved again with the periodic 
boundary condition. 



The observable quantities of clusters which we consider in this letter are the X-ray surface 
brightness and the spectral intensities due to the thermal and kinematic SZ effect at mm and 
submm bands. The X-ray surface brightness of clusters between the frequency bands v\ and V2 
for clusters located at z is given by the following integral along the line-of-sight (e.g., Rybicki, 
Lightman 1979): 



where a[vi,V2,T e ] is the X-ray emissivity in the corresponding frequency band, T e is the 
temperature of electron gas and n e is the number density of electrons. In practice we adopt the 



3. Fluxes in the X-ray, mm and submm bands 




(1) 
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Masai model (Masai 1984) for a[u±, V2>T e ] and thus include both metal line emissions and the 
thermal bremsstrahlung; the former becomes important for low temperature clusters. 

The spectral intensity due to the thermal SZ effect is written as (e.g., Rephaeli 1995): 
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where io = 2(/cb7cmb) 3 /(^pc) 2 , y is the Compton y-parameter : 



and g(x) is a function of x = h p v / '/cb^cmb: 
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In the above expressions, &b, n^e 0T, c, h p and Tcmb denote the Boltzmann constant, the electron 
mass, the Thomson cross section, the velocity of light, the Planck constant, and the temperature 
of the CMB, respectively. 



The corresponding change of the CMB temperature due to the thermal SZ effect is 
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The peculiar velocity of a cluster along the line-of-sight, v r , relative to the CMB rest frame 
(defined to be positive for a receding cluster) produces a kinematic SZ effect (Sunyaev, Zel'dovich 
1980): 
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where r is an optical depth of the cluster, and 

h{x) " (e- - iy 

Note that the last equality of equation ([?]) assumes that the temperature profile of the cluster is 
isothermal, the validity of which we will discuss further in §4.2. 

From the nine simulated clusters in total, we trace their progenitors at z = 1.0. We choose two 
clusters, A and B, so as to represent rich and poor clusters in our simulated sample, respectively. 
Table [l] summarizes their properties at z = and z = 1.0, in which the X-ray emission weighted 
temperature is defined as 

T e (r)a[T e ]n1(r)r 2 dr 

. (9) 
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cluster quantities 


clust 


er A 


clust 


er B 


z = 


z = 1 


z = 


z = 1 


r vir [Mpc] 


2.0 


0.80 


0.80 


0.40 


M vir [M ] 


1.3 x 10 15 


5.5 x 10 14 


1.2 x 10 14 


6.1 x 10 13 


-/Vgas(< ^vir) 


2186 


819 


14501 


6957 


Nbm(< r v ir) 


2801 


1191 


17379 


8482 


Lx,2-io [10 44 erg/sec] 


11 


10 


0.85 


1.2 


T x [keV] 


7.1 


5.9 


2.1 


1.4 


o"id [km/sec] 


1300 


1100 


587 


545 


/^spec 


1.4 


1.3 


1.0 


1.3 


L box [fc^Mpc] 


100 


25 



Table 1: Physical quantities characterizing the simulated clusters A and B at z = 0.0 and z = 1.0: 
the virial radius in physical lengths (r v i r ), the total mass within the virial radius (M v ; r ), the number 
of gas and dark matter particle within the virial radius (A" gas and Adm), the X-ray luminosity in 2- 
10 keV band (Lx,2— io), the X-ray emission weighted temperature (Tx), the 1-dimensional velocity 
dispersion of dark matter (a id), /3-parameter defined by f3 spec = /um p o^ D /fcBTx, anc ^ -^box is the 
comoving size of the simulation box. 

Figure [j] shows the projected profiles of the X-ray surface brightness, temperature decrement, and 
submm surface brightness at z = 0.05 and z = 1.0 for the two clusters, and dash-dotted lines show 
the /3-model profiles adopting their fitted parameters. The energy bands adopted in Figure |l| 
correspond to those of ASCA (2-10keV) in X-ray, NOBA (NOBeyama Bolometer Array; 150GHz) 
in mm, and SCUBA (350GHz) in submm. Note that we are mainly interested in quantifying 
the fractional uncertainties in the measurements of H$ and the peculiar velocities. In this sense, 
the absolute magnitude of the value of the y-parameter is not important; since our simulations 
did not include radiative cooling and other scale-dependent physical processes, the results are 
almost scalable and it seems unlikely that our results presented here change drastically for much 
larger clusters which have the y-parameter large enough to be observable in practice. We plan to 
come back to the issue with much larger simulations including various scale-dependent physical 
processes in due course. 

If the cluster gas strictly obeys the isothermal ^3— model: 

n e (r)=n e0 [l + (r/r c ) 2 ]- 3 ^ 2 , (10) 
then X-ray surface brightness (eq.|ij) and y-parameter (eq.||) should be given by 

(11) 
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and 

^ r(3/?/2 - 1/2) (k B T c /m c c 2 )a T n e or c 

y[ ) V7r r(3/3/2) [i + ^a^^Ac) 2 ] 3 ^ 2 - 1 / 2 ' 

where r c is the core radius of the cluster (in physical lengths), d\{z) is the angular diameter 
distance to the redshift z (e.g., Birkinshaw, Hughes, Arnaud 1991; Kobayashi, Sasaki, Suto 1996). 
Therefore we attempt the following fits to the projected profiles of clusters A and B: 

S x (9) oc [1 + (^ c ,x) 2 r %+1/2 , y{9) oc [1 + (0/e c , y ) 2 }-^ /2+1/2 - (13) 

The results of the fits to equations ( |io|) and (13) are summarized in Figure and also plotted 
in Figure [l] (thick solid and thin dashed lines for cluster A and B, respectively). The fits are 
performed in the range of £q < r < 2r v ; r for equation (|l0|), and of eo/d^ < 9 < 2r v i r /d\ for 
equation (13). Note that the core radii of both clusters decrease with redshift indicating that the 



clusters are still contracting from z = 1 to z = and thus are not in equilibrium at z = 1. 

While Figure |l] clearly illustrates that the separate fits to the /3-model predictions (eqs.|l 
and [13 1) are quite successful, the best-fitted values of r c (or 9 C ) and f3 change substantially 
depending on whether one uses n e (r), y(9), or Sx{9) for the fit. This should be ascribed to 
non-isothermality, asphericity and local dumpiness of clusters. A closer look at Figure § implies 
that the former is more important for clusters at z = (because the result is almost insensitive to 
the line-of-sight directions) while asphericity and local dumpiness dominate for clusters at z = 1. 

To understand the above, we plot the three dimensional profiles both for T e (r) and /0 gas (r) in 
Figure ||. The integrand of y is proportional to T e (eq. || ) while Sx is more weakly dependent on 
T e . Recall that the bolometric emissivity of the thermal bremsstrahlung is proportional to 
and the inclusion of line emission makes the dependence on T e even weaker for lower temperature 
clusters. Thus the gas temperature drop in the central regions, as shown in Figure ||, should 
increase both r c and j3 fitted to y(6), and to a lesser extent to Sx(9) compared with those fitted 
to n e (r), which is consistent with the systematic trends for clusters at z = 0. In Figure || different 
line-of-sight projections yield a large scatter of r c and (5 for clusters at z = 1 indicating that 
asphericity and dumpiness are more appreciable. 

Figure § displays the contours of |A/ mm | (eq.|§ at 150 GHz), AI su b mm (eq.|| at 350 GHz), 
Sx (eq. Q at 2 - 10 keV band), and the X-ray emission weighted temperature for clusters A and 
B. Apparently the contours of |A/ mm | and A/ su b mm are more extended than that of Sx, reflecting 
that the formers are essentially line-of-sight integrals of n e (r) rather than of n 2 (r) for Sx- Also 
they clearly exhibit that A/ mm and A/ su b m m do not suffer from (1 + z) -4 diminishing effect, which 
is quite suitable for cluster surveys at high redshifts compared with those in optical and X-rays. 



4. Estimating the Hubble Constant and Peculiar Velocity of Clusters 
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4.1. Hubble Constant 

The determination of the Hubble constant from cluster observation relies on the prediction of 
core radius of clusters of galaxies from the /3-model fit combining the SZ and X-ray fluxes (eqs.pfl, 



P and 



[Ar(0)/T]g bs r(3/3 fit , x - l/2)r(3/? fit , x /2) 2 
c ' X ' est S x (0) obs r(3^ fit)X )r(3/3 fiti x/2 - 1/2)2 

v ml^a\v x (\ + z),u 2 (l + z),T c ] 2 

4^/2(1 + ^4^2^ W> ^ 

(Silk, White 1978; Birkinshaw, Hughes, Arnaud 1991). 

The reliability of the estimated Hubble constant crucially depends on the relevance of the 
isothermal /3— model. Following Inagaki, Suginohara & Suto (1995), we examine the distribution 
of the parameter: 

f _ fc,X,6t Ha,est / 1E \ 

JH,ZD = 77 s - 77 \ lt> ) 

r c ,X,est{PXfit) -H0,true 

for all of the simulated clusters from three orthogonal directions at z = 0.0 and z = 1.0, where 
?"c,x,fit and /?x,fit are computed by fitting to the profile of Sx(6) (eq.fi"l|). The second equality 
in equation (|l5|) holds if the isothermal /3-model is exact and the angular diameter distance is 
approximated as ~ cz/Hq. For comparison, we also compute the similar quantity: 

p ^"c,n e ,fit /-, ^\ 

JH,SD = 77 r i (lOJ 

r c , n c , es t V Pn c , fit J 



using r C)Ilei fit and /3 nci fit fitted to the three-dimensional (spherically averaged) gas profile (|1CJ|). Note 
that fn,2Di rather than Jh,3Di should be considered as a measure of the observational uncertainties 
of Hq since only r C)X ,fit and Pxfit are directly estimated from the X-ray observations. 

Panels (a) to (d) in Figure || show the histograms of /h,2D and /h,3D for the nine clusters 
at z = 0.05 and z = 1.0. It is interesting to note that there is no significant systematic bias for 
Ih,2D at z = even though the values of r Cj x,fit and /3 x ,fit are fairly different from r c ricj fi t and 
Ai e , fit- In fact, the mean values of /#,2D and /h,3D are close to unity well within the la statistical 
errors. As first discussed by Inagaki, Suginohara, & Suto (1995), any temperature structure could 
bias the value of Hq estimated on the basis of equation (0). Although our simulated clusters 
do show some temperature structure, they do not seem to be strong enough to cause significant 
systematics. On the other hand, clusters at z = 1.0 exhibit the systematic underestimate bias of 
Hq by ~ 20%. This underestimate would result from the asphericity of the clusters because the 
values for the same clusters from three line-of-sight directions are significantly different. 



4.2. Peculiar Velocity 



The observed SZ flux AI U is contributed from both the thermal and kinematic SZ effects 
(eqs.||] and 0). By using the different spectral dependence of these effects, we can estimate the 
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peculiar velocity of the cluster (Sunyaev, Zel'dovich 1980; Rephaeli, Lahav 1991). Denote the total 
SZ flux in two bands, say v = v\ and i^, by 

AIj = i [y g(xj) - r(v r / 'c)h(xj)], (j = 1 and 2), (17) 

where Xj = h p Vj /A;b7cmb- While the y-parameter is observable, r is not unless the three- 
dimensional temperature profile is known independently. So one may estimate the peculiar velocity 
using the X-ray emission weighted temperature along the line-of-sight: 



T e (r)a[T e ]n e (r)dl 

(Tx) = J -^foo , (18) 

a[T e ]n1(r)dl 



as follows: 

v r ,est _ k B (T x ) ffQcjJAJg - g(x 2 )AI 1 ^ 
c m e c 2 h(xi)Al2 — h{x2)AI\ 

If clusters are isothermal, then y = t(/cb(7x)) / '(m e c 2 ) and v V)es t is identical to v T in equation (p^7[). 
For real clusters, however, the temperature profile is not strictly isothermal and estimates from 
equation (|l9| ) should be different from v r depending on the degree of non-isothermality. If we 
introduce a quantify the degree of non-isothermality of our simulated clusters: 

/, s I*52» (20) 

y m e c z 

then 

V r ,est = V r fv (21) 

So f v can be regarded as a correction factor which relates f rjes t to the correct peculiar velocity: 

We compute f v for all simulated clusters at z = 0.05 and z = 1 from three orthogonal 
line-of-sight directions (Figs. |||e and f). Compared with the histograms of /h,2D and /h,3D, that 
of f v is more centrally concentrated. This is because f v is sensitive only to the non-isothermal 
structure along the line-of-sight while largely free from either the existence of local dumpiness 
or asphericity unlike /h,2D and Jh,3D- Since (Tx) has more weights on the high-density central 
regions, /„ should be less than unity for clusters with central temperature drop as in our simulated 
ones. Panels (e) and (f) in Figure |5| indicates that this is the case. In summary, the estimates of 
v r are fairly reliable with overall fractional uncertainties ~ 10% at z ~ and ~ 20% at z = 1 



5. Conclusions and Discussion 

We have performed a series of numerical simulations with particular attention to the 
feasibility of estimating the Hubble constant H$ and the cluster peculiar velocity v r via multi-band 
observations. Let us briefly summarize our conclusions here. 
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(i) A conventional isothermal /3-model describes well our simulated clusters both in projected 
X-ray surface brightness and the SZ flux (Fig.|l]) as well as the three-dimensional gas density 
profile (Fig. |3|). This is consistent with the analytical (Makino, Sasaki, Suto 1998) and numerical 
work (Eke, Navarro, Frenk 1998) on the basis of the universal density profile of dark matter halo 
(Navarro, Frenk, White 1997). 

(ii) The best-fit values for the core radius and /^-parameter change significantly depending on 
whether one attempts to fit the X-ray surface brightness, the SZ flux or the gas density (Fig.||). 
This reflects the non-isothermal temperature profile, local dumpiness and aspherical gas structure. 

(iii) Provided that our simulated clusters constitute a representative sample of the observed 
clusters, the systematic and statistical errors in the estimates of Hq and v r are relatively small 
(Fig.[|); fractional uncertainties of the estimates of Hq amount to ~ 20% mainly due to the 
departure from the isothermal and spherical gas density distribution. Those of v r range ~ 10% at 
z ~ and ~ 20% at z = 1. Since the three-dimensional temperature structure is very difficult to 
reconstruct from the projected temperature map (compare Figs. ^ and |j), the correction for the 
non-isothermality is not easy. 

The similar temperature drop in the clusters discussed here was found earlier in simulations 
by Evrard (1990) and also reported in some observed clusters (e.g., Ikebe et al. 1997). Evrard 
(1990) ascribed this to the fact that earlier collapsed gas is only mildly shocked relative to the 
subsequent infalling gas. Although this is one reasonable interpretation, it could result simply 
from the shape of gravitational potential (Navarro, Frenk, White 1997; Makino, Sasaki, Suto 1998) 
or even from the lack of spatial resolution of simulations to represent the shock in the core regions. 
Therefore results from improved simulations are clearly needed; a larger number of gas and dark 
matter particles, ./V = 128 3 and even N = 256 3 , is necessary to resolve the cluster core reliably. 
Then it makes sense to incorporate effects of cooling and heating which have been neglected here 
mainly because they are not so important with the mass resolutions of our present simulations. 
We have considered only one fairly specific cosmological model, and it is interesting to examine 
the cosmological model dependence of the results. Nevertheless we hope that our current results 
have highlighted several potentially important implications for the future cluster observations in 
the X-ray, mm, and submm bands. 
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Fig. 1. — Projected profiles of X-ray surface brightness (top), the thermal SZ temperature 
decrement (middle), and the thermal SZ submm surface brightness (bottom) of clusters A and 
B at z = 0.05 (left) and z = 1.0 (right). The quoted error bars at each angular radius indicate the 
1-cr statistical dispersion computed from the 16 points on the projected circle with the radius. 
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Fig. 2. — The core radius (in physical lengths) and /3-parameter for clusters A (triangles) and B 
(circles) at z = (open symbols) and z = 1 (filled symbols). Different sizes of different symbols refer 
to the values from separate fits to n e (r), Sx and y-parameter (from three orthogonal line-of-sight 
directions for the latter two). 
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Fig. 3. — Spherically averaged profiles of gas temperature (upper) and gas density (lower) of clusters 
A and B at z = 0.05 (left) and z = 1.0 (right). The radius r is in physical lengths. 
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Fig. 4. — Projected views of clusters A (upper panels) and B (lower panels) at z = 1 and z ~ 0. A 
box of (4Mpc) 3 (in physical lengths) located at the center of each cluster is extracted. The X-ray 
emission-weighted temperature (Tx), X-ray surface brightness (Sx), and the SZ surface brightness 
at mm and submm bands (|A/ mm | and A/ su bmm) are plotted on the projected X-Y plane by 
integrating over the line-of-sight direction (Z). The X and Y coordinates are in the physical lengths 
at the corresponding redshift, and related to the angular coordinate 9 from the cluster center as 
9<Ia(z). At z ~ 0, oIa(z) can be replaced by the real distance to the cluster from the observer. 
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Fig. 5. — Distribution of /h,2D (left), /h,3D (middle), and f v (right) for all the simulated clusters 
(nine in total) at z = 0.05 and z = 1.0 viewed from three different line-of-sight directions. Different 
patterns of the histogram correspond to different clusters. The mean and la statistical errors are 
quoted in each panel. 
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